Superfluidity of metastable bulk glass para-hydrogen at low temperature 
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Molecular para-hydrogen (P-H2) has been proposed theoretically as a possible candidate for su- 
perfluidity, but the eventual superfluid transition is hindered by its crystallization. In this work, 
we study a metastable non crystalline phase of bulk P-H2 by means of the Path Integral Monte 
Carlo method in order to investigate at which temperature this system can support superfluidity. 
By choosing accurately the initial configuration and using a non commensurate simulation box, we 
have been able to frustrate the formation of the crystal in the simulated system and to calculate the 
temperature dependence of the one-body density matrix and of the superfluid fraction. We observe 
a transition to a superfluid phase at temperatures around 1 K. The limit of zero temperature is also 
studied using the diffusion Monte Carlo method. Results for the energy, condensate fraction, and 
structure of the metastable liquid phase at T = are reported and compared with the ones obtained 
for the stable solid phase. 



PACS numbers: 67.63.Cd, 67.80.5, 67.10.Ba, 02.70.Ss 
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I. INTRODUCTION 

Superfluidity and Bose-Einstein condensation (BEC) 
have been stunningly shown in metastable dilute alkali 
gases, magnetically confined at ultralow temperatureSf^ 
The extreme diluteness of these gases allows for the 
achievement of BEC with an almost full occupation of 
the zero-momentum state that has been possible to ob- 
serve and measure quite easily. This contrasts with the 
difficulties encountered in the measure of the condensate 
fraction in liquid *He, which amounts to only 8% at the 
equilibrium densityi^ However, liquid "^He is a stable su- 
perfluid below the lambda transition Tx = 2.17 K and 
is therefore a more easily accessible system. Before the 
blowup produced in the field of quantum fluids by the 
first experimental realization of BEC gases, liquid helium 
was the only paradigm of a superfiuid. From long time 
ago, there has been great interest in the search of super- 
fiuid condensed phases other than liquid helium. Spin- 
polarized atomic deuterium and tritium are predicted 
to be fermionic and bosonic liquids, respectively, in the 
limit of zero temperature.'^''* However, their experimen- 
tal study has proven to be very elusive due to their high 
recombination rate, and only the case of atomic hydro- 
gen, whose ground state is a gas, has been experimentally 
driven to its BEC statc;^ The next candidate for super- 
fluidity is molecular hydrogen, which has been studied for 
a long time<^ This seems a priori an optimal system due 
to its very light mass but it crystallizes at relatively high 
temperature as a consequence of the intensity of its inter- 
molecular attraction, without exhibiting any superfluid 
transition in the liquid phase. In the present work, we 
study the properties of metastable liquid or glass molec- 
ular hydrogen at very low temperatures using quantum 
Monte Carlo methods. 

In 1972, Ginzburg and Sobyanin^ proposed that any 
Bose liquid should be superfluid below a certain temper- 
ature T\, unless it solidifies at temperature Tf higher 



than T\. To give a first estimation of T\, they used the 
ideal Bose gas theory, obtaining 



Tx = 3.31 
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(1) 



where m is the atomic mass, g is the spin degeneracy, 
ks is the Boltzmann constant and p is the density of 
the system. Ginzburg and Sobyanin proposed molecular 
paro- hydrogen (P-H2) as a plausible candidate for super- 
fiuidity: being a spinless boson (g — 1) with a small 
mass, P-H2 should undergo a superfiuid transition at a 
relatively high temperature (according to Eq. ([1]) , Ta ~ 6 

K). _ 

The estimation of Ta, given by Eq. ([T]), is clearly in- 
accurate in the case of dense liquids because it cannot 
account for the observed dependence of Ta with the den- 
sity. In fact, Ta slightly decreases in liquid ^He when p 
increases, a manifestly opposite behavior to the increase 
with p^/^ given by the ideal gas formula ([T]). In order to 
provide a more reasonable estimation of Ta, Apenko^ pro- 
posed a phenomenological prescription for the superfluid 
transition, similar to the Lindemann criterion for classi- 
cal crystal melting. In this way, he was able to take into 
account quantum decoherence effects due to the strong 
interatomic potential and to relate the critical tempera- 
ture for superfluidity with the mean kinetic energy per 
particle above the transition. For P-H2, he concluded 
that Ta should vary between 1.1 K and 2.1 K, depending 
on the density of the system. 

Superfluid P-H2 is not observed in a stable form be- 
cause it crystallizes at temperature Tf = 13.8 K, which 
is signiflcantly higher than the expected Ta . Several stud- 
ies about crystal nucleation in P-H2 have been performed 
in order to understand if the liquid can enter a super- 
cooled phase, i.e., a metastable phase in which the liquid 
is cooled below its freezing temperature without forming 
a crystal. Maris et al^ calculated the rate r(T) of homo- 
geneous nucleation of the solid phase from the liquid as a 
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function of the temperature T, showing a maximum of T 
around T = 7 K and a rapid decrease at lower tempera- 
ture. This suggests that, if it would be possible to super- 
cool the liquid through the range where V is large, one 
might be able to reach a low-temperature region where 
the liquid is essentially stable. However, recent experi- 
ments have indicated that, even at T ~ 9 K, the rate 
of crystal growth is so high that the liquid phase freezes 
quickly into a metastable polymorph crystal. 

Even though several supercooling techniques have been 
proposed to create a metastable liquid phase in bulk p- 
H2r^"— none of them has proven so far to be success- 
ful and no direct evidence of superfluidity has been de- 
tected. However, there are evidences of superfluidity in 
several spectroscopic studies of small doped P-H2 clus- 
ters. In 2000, Grebenev et ali^ analyzed the rotational 
spectra of a linear carbonyl sulfide (OCS) molecule sur- 
rounded by 14 to 16 P-H2 molecules absorbed in a larger 
helium droplet, which fixes the temperature of the clus- 
ter. When P-H2 is immersed in a "^He droplet (T = 0.38 
K), the measured spectra shows a peak indicating the 
excitation of angular momentum around the OCS axis. 
On the other hand, if the small P-H2 cluster is put inside 
a colder ''He-^He droplet (T = 0.15 K), the peak dis- 
appears: the OCS molecule is then able to rotate freely 
inside the hydrogen cluster, pointing to the superfluid- 
ity of the surrounding P-H2 molecules. These results 
have been confirmed in a later experiment on small p- 
H2 clusters doped with carbon dioxide. From a precise 
analysis of the rotational spectra, it has been possible to 
measure the effective momentum of inertia of these small 
systems, and thus of their superfluid fraction ps, pro- 
viding a clear evidence of superfluidity in clusters made 
up of iV < 18 P-H2 molecules. These clusters are too 
small to extract reliable predictions of a metastable liq- 
uid phase and larger clusters would be desirable. To this 
end, Kuyanov-Prozument and Vileso\J^ have been able 
to stabilize liquid clusters with an average size of iV « lO''^ 
P-H2 molecules down to temperature T = 2 K, but they 
do not see any evidence of superfluidity. Other attempts 
of producing liquid P-H2 well below Tf (T = 1.3 K) are 
based on the generation of continuous hydrogen filaments 
of macroscopic dimensions 

The search for a superfluid P-H2 phase has been in- 
tense also from the theoretical point of view. The rather 
simple radial form of the p-H2-p-H2intcraction and the 
microscopic accuracy achieved by quantum Monte Carlo 
methods have stimulated a long-standing effort for de- 
vising possible scenarios where supercooled P-H2 could 
be studied. In practically all the cases, the search is fo- 
cused on systems of reduced dimensionality or in finite 
systems. Path integral Monte Carlo (PIMC) simulations 
of P-H2 films adsorbed on a surface with impurities ob- 
served superfluidity for some arrangements of these im- 
purities,— but these results were posteriorly questioned 
by other PIMC studiesii^ In a one-dimensional channel, 
like the one provided experimentally by narrow carbon 
nanotubes, it has been predicted a stable liquid phase 



in the limit of zero temperature. "'^^ The largest number of 
theoretical works have been devoted to the study of small 
clusters, both pure^S"— and doped with impurities 1^"— 
All these simulations show that P-H2 becomes superfluid 
below a certain temperature T — 1-2 K and that the su- 
perfluid fraction depends on the number of molecules of 
the cluster. When the cluster becomes larger than a cer- 
tain molecular number {N > 18-25), solid-like structures 
are observed and the superfluidity vanishes. 

In the present work, we deliver a PIMC study of a 
metastable glass/liquid phase at very low temperature. 
Our main purpose has been to determine for the flrst 
time at which temperature this metastable phase be- 
comes superfluid and the value of the superfluid density 
and condensate fraction close to this temperature. The 
simulations are carried out following schedules which are 
similar to the ones used in a recent study of a glass ^He 
phase evolving from a normal to a superfluid state (su- 
perglass).-'^'* Our results show that this transition tem- 
perature is T ~ 1 K, a value that is close to the Apenko 
estimatiorii^. and also close to the values observed in sim- 
ulations of small clusters. As a complementary aspect, 
we address the calculation of the equation of state of the 
metastable liquid P-H2 phase in the limit of zero temper- 
ature using the diffusion Monte Carlo (DMC) method. 
The simulation of the liquid phase in this limit is easier 
than at flnite temperature and therefore DMC is able to 
provide accurate information on its main energetic and 
structure properties. 

The rest of the paper is organized as follows. In Sec. II, 
we introduce the quantum Monte Carlo methods used in 
the study, the DMC and PIMC methods, and report spe- 
cific details on how the simulations are carried out. Sec. 
Ill contains the results of the equation of state, structure 
properties, and condensate fraction of metastable liquid 
P-H2 at zero temperature. PIMC results at finite tem- 
perature are reported in Sec. IV, and finally the main 
conclusions of the present work are discussed in Sec. V. 



II. QUANTUM MONTE CARLO METHODS 

The H2 molecule, which is composed of two hydrogen 
atoms linked by a covalent bond, is spherically symmet- 
ric in the para-hydrogen state (total angular momentum 
zero). The energy scale involved in electronic excitations 
(~ 10^ K) is orders of magnitude larger than the inter- 
molecular one 10^ K), thus modelling the P-H2-P-H2 
interaction by means of a radial pair-potential and con- 
sidering the molecules as point-like turns out to be jus- 
tified upon the condition of low or moderate pressures. 
In this work, we have chosen the well-known and com- 
monly used semiempirical Silvera-Goldman pair poten- 
tial. '^'^ This potential has proved to be accurate at low 
temperature and in the pressure regimes in which we are 
interested. 

The study in the limit of zero temperature has been 
performed with the DMC method. DMC is a first- 
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principles method which can access exactly the ground 
state of bosonic systems. It is a form of Green's Func- 
tion Monte Carlo which samples the projection of the 
ground state from the initial configuration with the op- 
erator exp [—{H — Eo)t]. Here, T-L is the Hamiltonian 



n 



2m 



N 



N 



i=i<j 



(2) 



Eq is a norm-preserving adjustable constant and t is the 
imaginary time. The simulation is performed by advanc- 
ing in r via a combination of diffusion, drift and branch- 
ing steps on walkers R (sets of 3A^ coordinates) repre- 
senting the wavefunction of the systemj^ The imaginary- 
time evolution of the walkers is "guided" during the drift 
stage by a guiding wavefunction (j)G, which is usually a 
good guess for the wavefunction of the system. This func- 
tion contains basic ingredients of the system as its sym- 
metry, phase and expected behaviors at short and long 
distances according to its Hamiltonian. Technically, (jjc 
allows importance sampling and thus reduces the vari- 
ance of the ground-state estimations. It is straightfor- 
ward to show that for the Hamiltonian H and any opera- 
tor commuting with it, the expectation value is computed 
exactly within statistical error. Other diagonal operators 
which do not fulfill this condition require of a special 
treatment, known as pure estimation, ■^^ which leads also 
for this case to unbiased results. 

The phase of the system is imposed within the typi- 
cal imaginary-time length by the guiding wave function. 
This property of the DMC method is here a key point if 
we are pursuing a investigation of the properties of the 
metastable liquid P-H2 phase. Then, for the liquid phase 
(j)G is taken in a Jastrow form 



N 



MR)= n 



(3) 



i=i<j 

with a two-body correlation functional 



f{r) = exp 



2 [r 



2 exp 



A 



(4) 



In order to compare the results obtained for the liquid 
phase with the ones corresponding to the stable hep solid 
we have carried out some simulations with a guiding wave 
function of Nosanow- Jastrow type 



N N 



(5) 



the set {rj} being the lattice points of a perfect hep lat- 
tice. Optimal values for the parameters entering Eq. (|4]) 
are b = 3.68 A, L = 0.2, A = 5.24 A, and A = 0.89 A 
for the liquid phase, and b = 3.45 A, L = 0.2, A = 5.49 
A, and A = 2.81 A for the solid one. The Nosanow term 



is chosen in Gaussian form, g{r) = exp(— 7r^). The den- 
sity dependence of the parameters in the Jastrow term 
is small, and neglected in practice when used in DMC, 
whereas the Nosanow term parameter 7 is optimized for 
the whole range of densities. We have used 256 and 180 
particles per simulation box for the liquid and hep solid 
phases, respectively. The number of walkers and time- 
step have been adjusted to reduce any bias coming from 
them to the level of the statistical noise. 

At finite temperature T, the microscopic description 
of the quantum system is made in terms of the ther- 
mal density matrix, pn{R' , R; (3) = (i?'|e^''^|i2), with 
/3 — {kBT)~^. The partition function Z, which allows 
for a full description of the properties of a given system, 
satisfies the relation 



Z = Tr(e 



» M 

/ J^dRipN (Ri,Ri+i;e) 1 



(6) 



that relies on the convolution property of the density 
matrix. In Eq. ([5]), e = P/M and the boundary condi- 
tion Rm+1 = Ri applies. The remarkable feature of Eq. 
(H]), on which PIMC is based, is that one can access to 
information at a temperature T by convoluting density 
matrices at higher temperature MT.^^ 

PIMC describes the quantum A^-body system consid- 
ering M different configurations Rj of the same system, 
whose sequence constitutes a path in imaginary time. 
This means that the A^-body quantum system is mapped 
onto a classical system of A^ ring polymers, each one com- 
posed by AI beads. The different beads can be thought as 
a way to describe the delocalization of the quantum parti- 
cle due to its zero-point motion. For sufficiently large M, 
one recovers the high-temperature density matrix, where 
it is legitimate to separate the kinetic contribution from 
the potential one (primitive action). In this way, it is 
possible by applying Eq. (jS]) to reduce the systematic er- 
ror due to the analytical approximation for below the 
statistical uncertainty. However, the primitive action is 
too simple to study extreme quantum matter and a bet- 
ter choice for the action is fundamental to reduce both 
the complexity of the calculation and ergodicity issues. 
To this end, we have used a high-order Chin action^i^ 
to obtain an accurate estimation of the relevant physi- 
cal quantities with reasonable numeric effort even in the 
low temperature regime, where the simulation becomes 
harder due to the large zero-point motion of particles. We 
have analyzed the dependence of the P-H2 energy on the 
parameter e and determined an optimal value e = 1/60 
K^^ for which the bias coming from the use of a finite e 
value is smaller than the characteristic statistical noise. 

A relevant issue one has to deal with when approach- 
ing the low temperature limit with PIMC simulations 
arises from the indistinguishable nature of the particles. 
In the path integral formalism, the exchanges between 
L different particles are represented by long ring poly- 
mers composed by L x M beads. If we study a bosonic 
system like P-H2 , the indistinguishability of the particles 
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does not affect the positivity of the integrand function in 
Eq. ^ and the symmetry of Z can be recovered by the 
samphng of permutations between the ring polymers. In 
the present study, we have used the Worm Algorithm^ 
which provides a very efficient sampling in permutation 
space. 

The key aspect of the worm algorithm is to work in 
an extended configuration space, containing not only the 
usual diagonal configurations made up of ring polymers, 
but also off-diagonal configurations which are character- 
ized by the presence of an open polymer (defined as the 
worm). By working with off-diagonal configurations, it is 
possible to sample the bosonic permutations by means of 
single-particle movements, like the swap update, whose 
acceptance rate can be made comparable to that of the 
other updates in the sampling of polymers. In order to 
fulfill this condition, it is important to optimize two pa- 
rameters of the worm algorithm. The first of them is C, 
which regulates the acceptance probability of the move- 
ments switching from diagonal to off-diagonal configura- 
tion and vice versa. In our simulations, we choose C to 
get the number of sampled off-diagonal configurations to 
about 65-70% of total number of configurations. In this 
way, the system is allowed to spend enough time both in 
off-diagonal configurations, where the sampling of per- 
mutations is done, and in diagonal configurations, where 
relevant observables such as the energy or the superfluid 
density are evaluated. The second one is Mg, which is 
the number of beads rebuilt in the swap update. The pa- 
rameter Ms has to be chosen as a compromise between 
a small value, which would make difficult the search of 
the partner of the worm in the swap movement, and a 
large value which would make difficult the reconstruction 
of the polymer once the partner has been chosen. In our 
simulations, we use Ms which maximizes the acceptance 
rate of the swap update: the typical value of Ms is about 
10% of the number of beads M. 



III. ZERO-TEMPERATURE RESULTS 

We have calculated the main properties of the 
metastable liquid and stable hep solid phases of P-H2. 
Our main goal has been to know the properties of a hy- 
pothetical bulk liquid phase and compare them with the 
ones of the stable solid. In order to achieve reliable esti- 
mations of liquid P-II2 it is crucial to work with a guiding 
wave function of liquid type, as we have discussed in the 
preceding Section. Within the typical imaginary-time 
length of our simulations we have not seen the forma- 
tion of any crystal structure, i.e., no signatures of Bragg 
peaks in the structure function S{k) have been registered 
so far. 

In Fig. [U we plot the DMC energies per particle of 
metastable liquid P-H2 as a function of the density. For 
comparison, we also report the results obtained for the 
hep crystal phase. Our hep energies are in close agree- 
ment with the ones reported in Ref. \^ using the same 
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FIG. 1. DMC energies per particle as a function of the density. 
Squares and circles correspond to the liquid and solid phases, 
respectively. Solid and dashed lines are the polynomial fits 
to the DMC energies for the liquid and solid, respectively. 
The diamond is the experimental energy of hep molecular 
hydrogen from Ref. \^ 

Silver a- Goldman potential. In the figure, we also show 
the experimental estimation at T = K from Ref. H^, 
E/N = —89.9 K, that lies a bit below of our results. 
This is again in agreement with previous DMC results^ 
which show that the experimental energy is, in absolute 
value, underestimated and overestimated by the Silver a- 
Goldman and Buck potential, respectively. Our results 
for both phases are well reproduced by the polynomial 
law 

(E/N)q and po being the equilibrium energy per parti- 
cle and equilibrium density, respectively. These equa- 
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FIG. 2. Pressure of the liquid (solid line) and solid (dashed 
line) P-H2 phases as a function of the density. Experimental 
points for the solid phaso^ are show as solid circles. 
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FIG. 3. Speed of sound of the liquid (solid line) and solid 
(dashed line) p-H2 phases as a function of the density. 



tions of state are shown in Fig. [T] with lines. The opti- 
mal parameters of the fits are: po — 0.026137(20) A^'^, 
{E/N)o = -87.702(37) K, ^ = 235(2) K, B = 140(10) 
K for the solid, and po = 0.023386(40) A^^^ {E/N)o = 
-76.465(51) K, A = 188(1) K, B ^ 131(10) K for the 
liquid. As expected, our DMC results shows that the 
solid phase is the stable one with a difference in energy 
per particle at the respective equilibrium points of ^ 10 
K, the equilibrium density of the liquid being 10 % 
smaller than the solid one. The same trend was observed 
in a DMC simulation of two-dimensional p-il2 , but there 
the differences were significantly smaller4^ It is worth 
noticing that about one half of the energy difference in 
the bulk systems comes from the decrease of the kinetic 
energy per particle going from the liquid to the solid: 
at density p = 0.03 A~^, it amounts to 93.3(1) and to 
89.5(1) K for the liquid and solid, respectively. 

From the equations of state (O, it is easy to deduce 
the pressure of the system at any density using the re- 
lation P{p) = p'^[d{E / N) I dp]. The results obtained for 
metastable liquid and stable solid phases are shown in 
Fig. [21 As one can see, at a given density the pressure 
of the liquid is larger than the one of the solid mainly 
because of the different location of the equilibrium densi- 
ties (P = 0). The results for the solid are compared with 
experimental data from Ref. |4^ showing a good agree- 
ment especially for not very large pressures. The density 
at which the function P{p) has a zero slope defines the 
spinodal point; beyond this limit the system is no more 
thermodynamically stable as a homogeneous phase. At 
this point, the speed of sound c{p) = [m~'^{dP/dp)Y^'^ 
becomes zero. Results for c(p) are shown for both phases 
in Fig. [31 The speed of sound decreases when the 
density is reduced and drops to zero at the spinodal 
point: (pc = 0.0176(1) A'^, = -12.6(5) MPa) and 
{pc = 0.0193(1) A-3^ Pc = -18.5(5) MPa) for liquid and 
solid, respectively. 

DMC produces also accurate results for the structure 
of the bulk system. In Fig. [31 we show results for the 
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FIG. 4. Top panel: Two-body radial distribution function 
of the liquid P-H2 phase at different densities: solid, long- 
dashed, short-dashed, and dotted lines stand for densities 
p = 0.0180, 0.0245, 0.0300, and 0.0340 A"^, respectively. Bot- 
tom panel: Static structure factor of the liquid phase. Same 
densities and notation than in the top panel. 

two-body radial distribution function g{r) of the liquid 
P-II2 phase for a set of densities. This function is pro- 
portional to the probability of finding two molecules sep- 
arated by a distance r. Increasing the density, the main 
peak becomes higher and moves to shorter interparticle 
distances; at least three peaks are observed. All these 
features point to the picture of a very dense liquid, with 
much more structure than in stable liquid "^He. In the 
same Fig. [4l we show results for the static structure 
factor S{k), related to g{r) by a Fourier transform. As 
one can see, the main peak increases quite fast with the 
density suggesting a highly structured metastable liquid. 
Nevertheless, we have not observed within the scale of 
the simulations the emergence of any Bragg peak which 
would point to formation of crystallites in the simulation 
box. In Fig. [5l we illustrate the comparison between S{k) 
for the liquid and solid systems at a density p = 0.0245 
A~'^, close to the equilibrium density of the liquid. The 
difference is the one expected between a liquid and a 
solid: oscillating function towards one at large k for the 
liquid and a sequence of Bragg peaks, corresponding to 
the hep lattice, for the sohd. 

One of the most relevant properties of a superfluid is 
the mean occupation of the zero-momentum state, i.e., 
the condensate fraction no. As it is well known, no can 
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be obtained from the asymptotic behavior of the one- 
body density matrix pi{r), 

no = lim pi{r) , (8) 

r— fcjo 

with pi (r) being obtained as the expectation vahie of the 
operator 

/ $(ri,...,ri+r, .■.,rjv) \ 
\ *(ri,...,r^) / ■ 

DMC resuhs for the condensate fraction of hquid P-H2 as 
a function of the density, obtained using the extrapolated 
estimator (there are no rehable pure estimators for non- 
diagonal operators), are shown in Fig. |6l The decrease 
of no with the density is well described by an exponen- 
tial decay (line in the figure). The strong interactions 
induced by the deep attractive potential well produce a 
big depletion of the condensate state. At the equilib- 
rium density, our estimation for the condensate fraction 
is no = 0.0037(7). This value is more than one order of 
magnitude smaller than the measured condensate frac- 
tion^ of liquid ''He at equilibrium (0.08). 

IV. SUPERFLUID TRANSITION 
TEMPERATURE 

One of the main goals of our work has been to de- 
termine the temperature at which a disordered phase of 
P-H2 becomes eventually superfluid. Recently, a similar 
approach has been used to study the superfluid properties 
of a glass phase of '*He,— a system that has been named 
superglass and that it has been argued to be related with 
some of the effects observed in torsional oscillator exper- 
iments on solid "^He. 

The first difficulty we have to deal with when inves- 
tigating computationally a disordered phase of P-H2 at 
low temperature is to provide a good equilibration of the 




FIG. 5. Static structure function of liquid and solid P-H2 at 
density p = 0.0245 A"^. The result for the hquid SL{k) (left 
scale) is shown with a solid line; the one for the hep solid 
Ss{k) (right scale) with a dashed line. 
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FIG. 6. Condensate fraction of metastable liquid p-B.2 as a 
function of the density. The points are the DMC results and 
the line is an exponential fit to them. 



system. The PIMC method, indeed, is aimed at study- 
ing the thermodynamic properties of quantum systems 
at thermal equilibrium. On the contrary, our purpose 
here is to study a configuration different from the one of 
minimum free energy, which for P-H2 at low temperature 
is the crystalline one. In order to do that, it is funda- 
mental to choose thoroughly the dimensions of the sim- 
ulation box and the number of particles, which must not 
be commensurate with any crystalline lattice. Also, it is 
important to choose a good initial configuration which 
evolves, as the Monte Carlo simulation goes on, towards 
a non crystalline phase which remains metastable for a 
number of Monte Carlo steps large enough to get a good 
statistics of the relevant quantities of the system. In this 
equilibration process, special attention must be paid to 
the thermalization of the polymers used within the PIMC 
formalism. A bad choice of initial conditions may cause 
the evolution of the system towards a configuration where 
the polymers representing each molecule are not allowed 
to spread and thus are not able to describe properly the 
zero-point motion of the molecules. This eventuality may 
represent a serious problem in our simulation, since we 
are mainly interested in the investigation of the super- 
fluid properties of P-H2. 

To check whether an equilibration scheme is efficient or 
not, it is important to monitor how the numerical estima- 
tions of the physical quantities change with the number 
of Monte Carlo steps. If we see that, as the simulation 
goes on, the computed variables do not show any evident 
trend but fluctuate around a certain value, we can con- 
clude that the system has reached the metastability. To 
check if this eventual metastable phase is crystalline or 
not, we can calculate the static structure factor S{k) and 
observe if it presents the Bragg peaks typical of a crystal 
configuration. 

We have used different technical schemes to get the 
desired metastable glass/liquid configuration. In many 
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FIG. 7. Static structure factor S{k) of metastable bulk P-H2 
at r = 10 K obtained with PIMC (solid line). This result 
is compared with the static structure factor of liquid P-H2 at 
zero temperature obtained with DMC (dotted line) 

cases, we were not able to stabilize this phase and after 
some time the liquid froze. Finally, we managed to de- 
vise a successful approach that is based in the following 
two steps. In the first part, we perform the simulation 
of a fictitious system of quantum particles with a mass 
equal to the one of the P-H2 molecule, but interacting 
through the ^He-'^He Aziz potential.^- Compared with 
the H2 Silver a- Goldman pair interaction, the Aziz poten- 
tial does not present a so deep attractive well and thus 
it is not able to freeze the system. Once this fictitious 
system is equilibrated, we change the Aziz interaction by 
the Silvera-Goldman one and equilibrate again towards 
the metastable P-H2 phase. At all the temperatures we 
consider in our study, we have verified that the super- 
fluid properties of P-H2 do not depend on the fact that 
permutations are allowed or not in the first step of the 
equilibration. To test this equilibration scheme, we have 
performed a simulation of iV = 100 P-H2 molecules inter- 
acting through the Silvera-Goldman potential^ inside a 
cubic box at the equilibrium density of the liquid phase 
at zero temperature (see Sec. Ill), p = 0.0234 A . For a 
preliminary test, we choose to perform the PIMC simu- 
lation at temperature T = 10 K, that is an intermediate 
temperature below the freezing temperature, where the 
liquid phase should be unstable, but above the estimated 
superfluid transition temperature, in order to make the 
simulation easier. It is worth noticing that after equili- 
bration the glass phase is better sampled, and thus crys- 
tallization is avoided, when the center of mass of all the 
polymers are moved simultaneously and accepted or not 
collectively by a single Metropolis step. 

In Fig. [71 we have plotted the static structure factor 
S(k) of the metastable phase and compared it with the 
same quantity computed for the liquid phase with DMC. 
We can see that the curve obtained at T = 10 K presents 
the first peak at the same k as the S{k) of the liquid and 
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FIG. 8. Static structure factor S{k) of metastable bulk P-H2 
at T = 2 K obtained with PIMC (solid line). For comparison, 
we also show the S{k) of liquid P-H2 at zero temperature 
obtained with DMC (dotted line). 

follows the same behavior up to the second maximum 
which is at fc ~ 4A . Even though the PIMC calcu- 
lation gives a peak which is higher and narrower than 
the peak obtained with DMC, and indicates that the 
PIMC configurations are slightly more structured than 
the DMC ones, we can conclude that our equilibration 
scheme is able to create a metastable liquid phase, at 
least in the range of intermediate temperature below the 
freezing point Tf and above the expected superfluid tran- 
sition T\. 

Since our main purpose is to localize the superfluid 
transition of this non crystalline phase, it is worth to test 
this equilibration scheme at lower temperatures, closer 
to the expected T\. For this reason, we have performed 
a simulation with = 90 P-II2 molecules at the same 
density, p — 0.0234 A , but at a lower temperature, 
T = 2 K. Once the mean value of the energy was stable, 
we computed the static structure factor S{k). The result 
is shown in Fig. |8l in comparison with the static structure 
factor of the zero-temperature liquid. As we can see, 
S{k) obtained in the PIMC simulation presents narrow 
maxima in the range of small k and is different from the 
typical S{k) of a liquid phase. However, these maxima 
tend to disappear at higher k and their height is much 
lower than the height of the Bragg peaks appearing in 
the <S'(fc) of a crystal. This indicates that the system 
simulated with PIMC has relaxed to a glass phase, which 
is structured at short range but lacks of the long-range 
coordination typical of the crystal structures. 

Even if a glassy configuration can make the diffusion 
of particles harder, the lack of long range coordination 
makes possible the appearance of off-diagonal long range 
order and it is worth to study the superfluid properties 
of this phase. To do that, we have studied the temper- 
ature dependence of the one-body density matrix pi{r). 
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FIG. 9. One-body density matrix pi(r) for glass P-H2 at 
density p — 0.0234 A and at different temperatures. T — 
0.7K (crosses), T = I. OK (open squares), T = I. OK and 
iV = 130 (triangles), T = 1.2K (diamonds), T = 1.5K (full 
squares), and T = 2.0 K (circles). 



The first simulation has been performed at T = 2 K, us- 
ing iV = 90 P-H2 molecules inside a cubic box at density 
p = 0.0234 A . The result for pi (r) obtained in this sim- 
ulation is shown in Fig. |9j at T = 2 K, we can clearly see 
an exponential decay of pi{r) at large r, indicating that 
Bose-Einstein condensation is not present in the system. 
In fact, we have noticed that in this calculation the swap 
update (the update responsible for bosonic exchanges in 
our PIMC sampling) has a very low acceptance rate and 
does not allow the formation of long-permutation cycles 
with a non-zero winding number. Nevertheless, one can 
think that the low acceptance of the swap update is a 
consequence of the difficulties in the sampling of the co- 
ordinates due to the strength of the intermolecular po- 
tential. We may therefore suspect that the system re- 
mains stuck in a configuration without permutation be- 
cause of sampling issues. To be sure about our result, 
we have performed another simulation of the same sys- 
tem but starting from an initial configuration presenting 
a non-zero winding number. To create this initial con- 
figuration, we have allowed particles to permute even in 
the fictitious simulation used to equilibrate the system. 
When we start the PIMC simulation of P-H2 from the 
permutated configuration, we see that the percentage of 
particles involved in bosonic exchanges tends to decrease 
and, at the end of the equilibration, the system has re- 
laxed to a phase presenting zero winding number, i.e., 
the superfluid density is zero. This last result confirms 
our conclusion that the P-II2 glass we simulate is not 
superfluid at T = 2 K. 

In Fig. [9l we have also shown pi (r) estimated at other 
temperatures. At each of the temperatures studied, we 
have performed simulations starting both from a per- 
mutated and a non permutated configuration, observing 



that once the system has been equilibrated the results for 
Pi(r) do not depend on the initial configuration. From 
the comparison of the curves at different temperatures, 
we can easily see a change of the behavior of pi at large 
r as the temperature decreases: this indicates that, at 
temperatures close to T = 1 K, the system presents a 
transition to a superfluid phase presenting off-diagonal 
long range order. The condensate fraction at low tem- 
perature is no ^ 3 X 10^*, appreciably smaller than the 
DMC estimation for the liquid phase at T = 0. The ob- 
servation of a flnite value for the condensate fraction at 
T — 1 K agrees with the measure of a flnite value for the 
superfluid density. Our results show that the superfluid 
density, derived from the winding number estimation, is 
zero within our numerical uncertainty for T > 1.2 K and 
at T = 1 K is already Ps/p = 0.36 ± 0.08. In order to 
give a more precise estimation of the superfluid transition 
temperature Tx , it would be necessary to perform a perti- 
nent flnite-size scaling study. However, the achievement 
of the metastable state is quite hard for systems made 
up of more than 100 molecules. In our study, we have 
been able to stabilize the amorphous phase for a system 
made up of 130 P-H2 molecules at T = 1 K. The result for 
pi (r) is also shown in Fig. [9] We notice that this result 
is in agreement with the one calculated for the smaller 
system of 90 p-il2 molecules at the same temperature, 
supporting our conclusion that the system has undergone 
a superfluid transition. At higher temperatures, instead, 
the system made up of 130 P-H2 molecules relaxes to a 
crystalline phase and it has not been possible to calculate 
pi (r) for the amorphous conflguration. This result seems 
to indicate that the disordered phase is somehow "more 
stable" at lower temperature, when the p-il2 molecules 
begin to permutate. A similar behavior has also been 
found in PIMC simulations of small P-II2 clustersi^ 



V. CONCLUSIONS 

We have carried out extensive quantum Monte Carlo 
calculations of P-II2 at temperatures well below its solid- 
iflcation point. Our interest has been to know better the 
properties of the metastable liquid/glass phase at very 
low temperatures and to determine where the superfluid 
transition is expected to appear. In the limit of zero tem- 
perature we have used the DMC method, which is a very 
efficient tool to sample this metastable phase through 
the use of a proper guiding wave function. The results 
point to a very structured liquid with a large depletion 
of the condensate fraction, significantly larger than in 
stable liquid "^He. 

Our estimation of Ta ~ 1 K is slightly smaller than the 
prediction obtained using a phenomenological approach^ 
for which, at the density p — 0.0234 A studied in our 
simulations, the transition temperature is estimated to be 
Ta ^ 1.7K. It is also interesting to notice that our result 
for T\ is quite close to the temperatures at which, accord- 
ing to PIMC simulations, superfiuid effects should appear 
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in small P-H2 clustersi^i^ These calculations show that 
clusters made of TV < 20 P-H2 molecules exhibit a non- 
zero superfluid fraction below T ^ 2 K. This transition 
temperature depends on the dimension of the cluster, de- 
creasing when the number of molecules increases. How- 
ever, it is difficult to make hypothesis on the superfluid 
behavior of large enough P-H2 systems from the simula- 
tion of small clusters, because the calculated superfluid 
fraction is significantly depressed when the number of 
molecules becomes iV > 30. This unexpected behavior of 
Ps with N has been explained by relating the changes in 
the superfiuid properties to structural changes that make 



the molecules arrange according to a solidlike configura- 
tion when the dimension of the cluster becomes large. 
In our simulation of bulk glass P-H2, we have been able 
to frustrate crystallization with an efficient equilibration 
of the system and to measure finite values of both the 
condensate fraction and superfiuid density. 
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